% function to detrend, de-bias, calibrate, fft, and cross & auto power calculation
% References Bendat & Piersol, 1971, Metronix 1997
% included sub functions parzen.m and doparzen.m in this file
% to save computing time 30.4.3
% checked, checked & checked the parzening pl do not do it again
% tried to bypass parzen averaging with simple median - very bad results 22.5.3
% Robust band average good results 5.6.3
% latest date 16.6.3
% ~10 years later
% changes made from the original settlin.m
% to get it work with obs time series
%latest date Feb 22, 2011
%#eml

function[SPMatrix,ProcDef] = GetSpectra_h(TS,ProcDef)

nchannel = TS.head.nchannel;
block = ProcDef.block;
TLine1 = ProcDef.TLine1;
TLRad1 = ProcDef.TLRad1;
TLFreq1 = ProcDef.TLFreq1;
nfrq = ProcDef.nfrq;
Array = ProcDef.Array;
[TIndex,len] = GetSelArray(ProcDef);
ProcDef.len=len;
ProcDef.SelStacks = ones([1,len]);
ProcDef.Weights = ones([nfrq,len]);
Channel1 = zeros([ProcDef.block,nchannel]);
Channel = zeros([ProcDef.block/2,nchannel]);

if len < 1,
   fprintf('At least one stack required - error in GetSpectra\n');
   return;
end;


%------------native declarations---------------
disp('GetSpectra -> Now making matrices of cross and auto spectra....');

WI = hanning(block);


SPMatrix = zeros([nfrq len nchannel nchannel]);

% fff = len/20;
% aaa = 1;
% fprintf('0                 100percent\n');
for i = 1:len,
   
%    if i/fff > aaa,
%       fprintf('#');
%       aaa = aaa+1;
%    end;
   
   
% FFT, calibration etc.


for u = 1:nchannel,
   K = CorTre(TS.matrix(TIndex(i,1):TIndex(i,2),u),0);
	Channel1(:,u) = fft(K.*WI)/block;
   %Channel(:,Array(u)) = Channel1(2:(block/2)+1,u).*conj(CalData(u,1:block/2)'); %
   Channel(:,Array(u)) = Channel1(2:(block/2)+1,u); %
end; 	
% OOOPS the ' transpose made all the probelm. Solved it with another conj -30,10,02 (b'day eve)
% auto and cross spectra---------------------------
for ii = 1:nfrq,
   par = parzen1(TLRad1(ii)); % Bit inefficient (as it is called many times to do the same job)
   SPMatrix(ii,i,:,:) = CompuSPM1(Channel,par,TLine1(ii),TLRad1(ii));
end;
%------------------------------------------------
%clear Channel;

end;
fprintf('\n Over..Switching to main\n');

%helper functions------------------------------

function[SPM] = CompuSPM1(Channel, par,TL,TR),

a = size(Channel);
nchannel = a(2);
SPM = zeros([nchannel,nchannel]);

for i = 1:nchannel,
   for j = 1:i,
%talkal = sum(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i))); % just summing 
talkal = sum(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i)).*par); % summing w.r.t parzen window%
%talkal = rbpz(Channel(TL-(TR-1):TL+(TR-1),j).*conj(Channel(TL-(TR-1):TL+(TR-1),i))); % Robust Band averaging
      if i ~= j,
       SPM(i,j) = real(talkal);
       SPM(j,i) = (-1)*imag(talkal); 
      elseif i == j,
       SPM(i,j) = talkal;
     end;
  end;
end;

function[PF] = parzen1(radius)
 
PF = zeros([1,radius]);

if radius > 1,

%---- compute norm-factor for parzening 

      nFactor = 1.0; 
		PF(1) = 1.0;
      for line = 2:radius,

        u = (line/(radius+1))*pi;
        v = (sin(u)/u)^2.0;
        PF(line) = v;
        nFactor = nFactor+(2*v*v);
      end;
      nFactor =sqrt(nFactor);

%      -- parzening line 1 

      PF(1) = 1.0/nFactor;

%     -- parzening line 2 to radius      
	for line = 2:radius,
        PF(line) = PF(line)/nFactor;
    end;
    ParInt = sum(PF(2:radius))*2 + PF(1);
    PF = PF/ParInt;
else
    PF = 1;
 end;
   PF= [PF(radius:-1:2) PF]';

   
   
   %----trying to make parzen average a robust est 25.5.3''
   
% function[C] = rbpz(S)
% 
% sig_m = 1.483*median(abs(S-median(S))); % median absolute deviation 
% c_m = 1.5*sig_m;								% Scale
% k = find(abs(S) <= c_m);					% assigning Huber weights 
% w_m(k) = 1;
% q = find(abs(S) > c_m);
% w_m(q) = c_m./abs(S(q));					
% C = sum(S.*w_m')/sum(w_m);  			% Weighted X spectra
